setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
library(openxlsx)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210316LIBD.eQTL处理/53K.add.GWAS.eQTL.DEG.motif.for.analysis.txt",header=T,sep="\t")
file=file[!duplicated(file$id),]	#对TF去重
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

asm2=read.table("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/220520ASMs_anno.hg19_multianno.csv",head=T,sep=",")	#220K
asm2$unitID=paste(asm2$Chr,asm2$Start,asm2$Ref,asm2$Alt,sep=":")
asm=read.csv("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/869727.all.snp.vaf.up.down.20210321.csv",head=T)
asm=asm[asm$unitID %in% asm2$unitID,]
asm=asm[asm$unitID %in% filea$unitID,]
asmdata=filea1[filea$unitID %in% asm$unitID,]

sci=read.csv("E:/0 公共数据库差异情况/2018 Science supp/aat8127_Table_S1a_DGE.csv",header=T)		###DEG数据合并
sci_scz=sci[sci$SCZ.fdr<0.1,]	#6324 sig gene
sci_bd=sci[sci$BD.fdr<0.1,]		#2010 sig gene
sci_scz=data.frame(symbol=sci_scz$gene_name,SCZ.log2FC=sci_scz$SCZ.log2FC,SCZ.2018sci.fdr=sci_scz$SCZ.fdr)
sci_bd=data.frame(symbol=sci_bd$gene_name,BD.log2FC=sci_bd$BD.log2FC,BD.2018sci.fdr=sci_bd$BD.fdr)
#sci_scz$group.2018sci.scz="SCZ"
#sci_bd$group.2018sci.BD="BD"

cmc.nosva=read.xlsx("E:/0 公共数据库差异情况/CMC/DEG_nosva.xlsx",sheet=1)
cmc.nosva=cmc.nosva[cmc.nosva$adj.P.Val<0.1,]	#1265 sig gene
cmc.nosva=cmc.nosva[!cmc.nosva$MAPPED_genes==".",]
cmc.nosva=data.frame(symbol=cmc.nosva$MAPPED_genes,cmc.logFC=cmc.nosva$logFC,cmc.nosva.fdr=cmc.nosva$adj.P.Val)
cmc.nosva$cmc.nosva.group="CMC.nosva"

#deg=merge(sci_scz,sci_bd,by="symbol",all=T)
#deg=merge(deg,cmc.nosva,by="symbol",all=T)

DEG=cmc.nosva
con=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201207/con_genotype.random3.anno.hg19_multianno.csv",header=T)
conid=unique(unlist(strsplit(con$Gene.refGene,";")))
conid=conid[!conid=="."]

c=length(intersect(as.character(DEG$symbol),as.character(conid)))
d=length(conid)-c

upDEG=DEG[DEG[,2]>0,]$symbol
downDEG=DEG[DEG[,2]<0,]$symbol
tmpstr=intersect(as.character(DEG$symbol),as.character(conid))

result=data.frame(matrix(,nrow = 8,ncol =10))
names(result)=c("Term","overlap.num","OR","P.value","up.DEG.num","OR.up","P.value.up","down.DEG.num","OR.down","P.value.down")	#分别记录总的overlap、upDEG、downDEG的富集
result[1,]=c("control.150K",c,NA,NA,length(intersect(tmpstr,upDEG)),NA,NA,length(intersect(tmpstr,downDEG)),NA,NA)

fileaid=unique(unlist(strsplit(filea$Gene.refGene,";")))
fileaid=fileaid[!fileaid=="."]
a=length(intersect(as.character(fileaid),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(fileaid))
tmp=fisher.test(matrix(c(a,c,length(fileaid)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(fileaid)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(fileaid)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[2,]=c("117012",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)

filea1id=unique(unlist(strsplit(filea1$Gene.refGene,";")))
filea1id=filea1id[!filea1id=="."]
a=length(intersect(as.character(filea1id),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(filea1id))
tmp=fisher.test(matrix(c(a,c,length(filea1id)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(filea1id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(filea1id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[3,]=c("61725",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)


fileid=unique(unlist(strsplit(file$Gene.refGene,";")))
fileid=fileid[!fileid=="."]
a=length(intersect(as.character(fileid),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(fileid))
tmp=fisher.test(matrix(c(a,c,length(fileid)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(fileid)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(fileid)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)

file1id=unique(unlist(strsplit(file1$Gene.refGene,";")))
file1id=file1id[!file1id=="."]
a=length(intersect(as.character(file1id),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(file1id))
tmp=fisher.test(matrix(c(a,c,length(file1id)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(file1id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(file1id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)

file2id=unique(unlist(strsplit(file2$Gene.refGene,";")))
file2id=file2id[!file2id=="."]
a=length(intersect(as.character(file2id),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(file2id))
tmp=fisher.test(matrix(c(a,c,length(file2id)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(file2id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(file2id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)

file3id=unique(unlist(strsplit(file3$Gene.refGene,";")))
file3id=file3id[!file3id=="."]
a=length(intersect(as.character(file3id),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(file3id))
tmp=fisher.test(matrix(c(a,c,length(file3id)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(file3id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(file3id)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)

asmid=unique(unlist(strsplit(asmdata$Gene.refGene,";")))
asmid=asmid[!asmid=="."]
a=length(intersect(as.character(asmid),as.character(DEG$symbol)))
tmpstr=intersect(as.character(DEG$symbol),as.character(asmid))
tmp=fisher.test(matrix(c(a,c,length(asmid)-a,d),nrow = 2,byrow = T))
tmp2=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,5]),length(asmid)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,5])),nrow = 2,byrow = T))
tmp3=fisher.test(matrix(c(length(intersect(tmpstr,upDEG)),as.numeric(result[1,8]),length(asmid)-length(intersect(tmpstr,upDEG)),length(conid)-as.numeric(result[1,8])),nrow = 2,byrow = T))
result[8,]=c("13649",a,tmp$estimate,tmp$p.value,length(intersect(tmpstr,upDEG)),tmp2$estimate,tmp2$p.value,length(intersect(tmpstr,downDEG)),tmp3$estimate,tmp3$p.value)

write.csv(result,"./20210323.DEG富集/CMC.enrich.con.150K.random3.csv")

###					以下对MAGMA方法分析GWAS数据的结果进行分析
BD.GWAS.MAGMA=read.table("E:/0 公共数据库差异情况/MAGMA.GWAS.疾病风险基因/BD_GWAS_MAGMA.result.txt",head=T,sep="\t")
#BD.GWAS.MAGMA$fdr=p.adjust(BD.GWAS.MAGMA$P,method = "BH")
BD.GWAS.MAGMA=BD.GWAS.MAGMA[BD.GWAS.MAGMA$P<0.05,]
BD.GWAS.MAGMAid=unique(as.character(BD.GWAS.MAGMA$V6))

BDSCZ.GWAS.MAGMA=read.table("E:/0 公共数据库差异情况/MAGMA.GWAS.疾病风险基因/BDSCZ_GWAS_MAGMA.result.txt",head=T,sep="\t")
BDSCZ.GWAS.MAGMA=BDSCZ.GWAS.MAGMA[BDSCZ.GWAS.MAGMA$P<0.05,]
BDSCZ.GWAS.MAGMAid=unique(as.character(BDSCZ.GWAS.MAGMA$V6))

SCZ.GWAS.MAGMA=read.table("E:/0 公共数据库差异情况/MAGMA.GWAS.疾病风险基因/CLOZUK_PGC2_GWAS_MAGMA分析结果.txt",head=T,sep="\t")#基于CLOZUK+PGC2的精分GWAS数据
SCZ.GWAS.MAGMA=SCZ.GWAS.MAGMA[SCZ.GWAS.MAGMA$P<0.05,]
SCZ.GWAS.MAGMAid=unique(as.character(SCZ.GWAS.MAGMA$V6))

symbols=BD.GWAS.MAGMAid

con=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201207/con_genotype.random1.anno.hg19_multianno.csv",header=T)
conid=unique(unlist(strsplit(con$Gene.refGene,";")))
conid=conid[!conid=="."]

c=length(intersect(as.character(symbols),as.character(conid)))
d=length(conid)-c

result=data.frame(matrix(,nrow = 7,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.150K",c,NA,NA)

fileaid=unique(unlist(strsplit(filea$Gene.refGene,";")))
fileaid=fileaid[!fileaid=="."]
a=length(intersect(as.character(fileaid),as.character(symbols)))
tmp=fisher.test(matrix(c(a,c,length(fileaid)-a,d),nrow = 2,byrow = T))
result[2,]=c("117012",a,tmp$estimate,tmp$p.value)

filea1id=unique(unlist(strsplit(filea1$Gene.refGene,";")))
filea1id=filea1id[!filea1id=="."]
a=length(intersect(as.character(filea1id),as.character(symbols)))
tmp=fisher.test(matrix(c(a,c,length(filea1id)-a,d),nrow = 2,byrow = T))
result[3,]=c("61725",a,tmp$estimate,tmp$p.value)

fileid=unique(unlist(strsplit(file$Gene.refGene,";")))
fileid=fileid[!fileid=="."]
a=length(intersect(as.character(fileid),as.character(symbols)))
tmp=fisher.test(matrix(c(a,c,length(fileid)-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

file1id=unique(unlist(strsplit(file1$Gene.refGene,";")))
file1id=file1id[!file1id=="."]
a=length(intersect(as.character(file1id),as.character(symbols)))
tmp=fisher.test(matrix(c(a,c,length(file1id)-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

file2id=unique(unlist(strsplit(file2$Gene.refGene,";")))
file2id=file2id[!file2id=="."]
a=length(intersect(as.character(file2id),as.character(symbols)))
tmp=fisher.test(matrix(c(a,c,length(file2id)-a,d),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value)

file3id=unique(unlist(strsplit(file3$Gene.refGene,";")))
file3id=file3id[!file3id=="."]
a=length(intersect(as.character(file3id),as.character(symbols)))
tmp=fisher.test(matrix(c(a,c,length(file3id)-a,d),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value)

write.csv(result,"./20210323.DEG富集/MAGMA.BD.enrich.con.150K.random1.csv")

library(ggplot2)
scz=read.csv("./20210317.DEG富集/SCZenrich.con.150K.random1.csv",head=T)
bd=read.csv("./20210317.DEG富集/BDenrich.con.150K.random1.csv",head=T)
cmc=read.csv("./20210317.DEG富集/CMCenrich.con.150K.random1.csv",head=T)

library(ggplot2)							###总体来看DEG的富集结果
scz=read.csv("./20210323.DEG富集/sci_scz.enrich.con.150K.random1.csv",head=T)[6,]
bd=read.csv("./20210323.DEG富集/sci_bd.enrich.con.150K.random1.csv",head=T)[6,]
cmc=read.csv("./20210323.DEG富集/CMC.enrich.con.150K.random1.csv",head=T)[6,]
MAGMA.scz=read.csv("./20210323.DEG富集/MAGMA.SCZ.enrich.con.150K.random1.csv",head=T)[6,]
MAGMA.sczbd=read.csv("./20210323.DEG富集/MAGMA.BDSCZ.enrich.con.150K.random1.csv",head=T)[6,]
MAGMA.bd=read.csv("./20210323.DEG富集/MAGMA.BD.enrich.con.150K.random1.csv",head=T)[6,]

scz=data.frame(counts=scz$overlap.num,P.value=scz$P.value,OR=scz$OR)
bd=data.frame(counts=bd$overlap.num,P.value=bd$P.value,OR=bd$OR)
cmc=data.frame(counts=cmc$overlap.num,P.value=cmc$P.value,OR=cmc$OR)
MAGMA.scz=data.frame(counts=MAGMA.scz$overlap.num,P.value=MAGMA.scz$P.value,OR=MAGMA.scz$OR)
MAGMA.sczbd=data.frame(counts=MAGMA.sczbd$overlap.num,P.value=MAGMA.sczbd$P.value,OR=MAGMA.sczbd$OR)
MAGMA.bd=data.frame(counts=MAGMA.bd$overlap.num,P.value=MAGMA.bd$P.value,OR=MAGMA.bd$OR)


library(ggrepel)
library(ggplot2)
scz$group="SCZ"
bd$group="BD"
cmc$group="CMC"
MAGMA.scz$group="MAGMA.scz"
MAGMA.sczbd$group="MAGMA.sczbd"
MAGMA.bd$group="MAGMA.bd"

deg=rbind(scz,bd,cmc,MAGMA.scz,MAGMA.sczbd,MAGMA.bd)
deg$group1=ifelse(deg$P.value<0.05,"S","N.S")

deg$log2Pvalue=-log2(deg$P.value)

ggplot(deg,aes(log2Pvalue,OR,group=group1,color=group1))+geom_point(aes(size=3))+
  scale_color_manual(values = c('#696969','#E64B35'))+
  theme_light(base_size = 15)+geom_vline(xintercept = -log2(0.05),color="red",linetype="dashed")+
    geom_text_repel(data = deg,aes(label = paste(group,counts,sep="\n")),size = 5)+
    theme(panel.grid.minor = element_blank())+
    xlab("-log2(P value)")#+ theme(legend.position = "none")

library(ggplot2)			###看差异表达的基因中上调性和下调性的富集情况
scz=read.csv("./20210323.DEG富集/sci_scz.enrich.con.150K.random1.csv",head=T)[6,]
bd=read.csv("./20210323.DEG富集/sci_bd.enrich.con.150K.random1.csv",head=T)[6,]
cmc=read.csv("./20210323.DEG富集/CMC.enrich.con.150K.random1.csv",head=T)[6,]

scz.up=data.frame(Term="SCZ",counts=scz$up.DEG.num,OR=scz$OR.up,P.value=scz$P.value.up,group="Up.DEG")
scz.down=data.frame(Term="SCZ",counts=scz$down.DEG.num,OR=scz$OR.down,P.value=scz$P.value.down,group="Down.DEG")
bd.up=data.frame(Term="BD",counts=bd$up.DEG.num,OR=bd$OR.up,P.value=bd$P.value.up,group="Up.DEG")
bd.down=data.frame(Term="BD",counts=bd$down.DEG.num,OR=bd$OR.down,P.value=bd$P.value.down,group="Down.DEG")
cmc.up=data.frame(Term="CMC",counts=cmc$up.DEG.num,OR=cmc$OR.up,P.value=cmc$P.value.up,group="Up.DEG")
cmc.down=data.frame(Term="CMC",counts=cmc$down.DEG.num,OR=cmc$OR.down,P.value=cmc$P.value.down,group="Down.DEG")

deg=rbind(scz.up,scz.down,bd.up,bd.down,cmc.up,cmc.down)
deg$group1=ifelse(deg$P.value<0.05,"S","N.S")
deg$term=paste(deg$Term,deg$group,sep=" ")
deg$log2Pvalue=-log2(deg$P.value)

ggplot(deg,aes(log2Pvalue,OR,group=group1,color=group1))+geom_point(aes(size=3))+
  scale_color_manual(values=c('#696969','red'))+
  theme_light(base_size = 15)+geom_vline(xintercept = -log2(0.05),color="red",linetype="dashed")+
  geom_text_repel(data = deg,aes(label = paste(term,counts,sep="\n")),size = 5)+
  theme(panel.grid.minor = element_blank())+
  xlab("-log2(P value)")#+ theme(legend.position = "none")

